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Abstract 

In this work we simulate a viscous hydrodynamical model of non-central Au-Au collisions in 2+1 
dimensions, assuming longitudinal boost invariance. The model fluid equations were proposed by 
Ottinger and Grmela Freezeout is signaled when the viscous corrections become large relative 
to the ideal terms. Then viscous corrections to the transverse momentum and differential elliptic 
flow spectra are calculated. When viscous corrections to the thermal distribution function are not 
included, the effects of viscosity on elliptic flow are modest. However, when these corrections are 
included, the elliptic flow is strongly modified at large pT- We also investigate the stability of the 
viscous results by comparing the non-ideal components of the stress tensor (vr*-^) and their influence 
on the V2 spectrum to the expectation of the Navier-Stokes equations {ir^^ = —r] (diUj)). We argue 
that when the stress tensor deviates from the Navier-Stokes form the dissipative corrections to 
spectra are too large for a hydrodynamic description to be reliable. For typical RHIC initial 
conditions this happens for rj/s ^ 0.3. 
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I. INTRODUCTION 



A. Motivation 



One of the first and most exciting observations from the Relativistic Heavy Ion Collider 
(RHIC) was the very strong elliptic flow in non-central collisions jH, [3|. The elliptic flow is 
quantified by the anisotropy of particle production with respect to the reaction plane f 2 

--(4^4). (1.1) 

and can be measured as a function of pt, rapidity, centrality and particle type. 

The adopted interpretation of the measurements is that the medium responds as a 
fluid to the differences in pressure gradients in the x and y directions. The fluid then 
expands preferentially in the reaction plane and establishes the observed momentum space 
anisotropy. This hydrodynamic interpretation is supported by the qualified success of ideal 
hydrodynamic models in describing a large variety of data over a range of colliding systems 
and energies jl, [H], 0, 0, Isj . Nevertheless, the hydrodynamic interpretation of the flow results 
is not unassailable. A back of the envelope estimate of viscous corrections to hydrodynamic 
results jsl suggests that viscous corrections are actually rather large, i.e., the mean free path 
is comparable to the system size These estimates are best conveyed in terms of the 

shear viscosity to entropy ratio, 77/5. The conditions for partial equilibrium at RHIC are so 
unfavorable that at unless 77/5 is small (say 0.5 or less), it is difficult to imagine that the 
medium would participate in a coordinated collective flow. 

From a theoretical perspective, it is difficult to reliably estimate t^/s in the vicinity of the 
QCD phase transition where the system is strongly coupled. Lattice QCD measurements of 
transport are hard (perhaps impossible 11, though recent efforts have lead to estimates 



which are not incompatible with the hydrodynamic interpretation of RHIC results |l3l . Il4 |. 
In a strict perturbative setting (where the quasi-particle picture is exact) 77/5 is large ~ l/f?^. 
Nevertheless an extrapolation of weak coupling results to moderate coupling also leads to 



an r^/s which is perhaps reconcilable with the hydrodynamic interpretation |15l. Il6|. Finally, 



these perturbative estimates should be contrasted with A/" = 4 Super Yang Mills at strong 
coupling, where r7/sisl/47r[l7,18]. Although A/" = 4 SYM is not QCD, the calculation was 
important because it showed that there is at least one theory where r^/s is sufficiently small 
that collective phenomena would be observed under conditions similar to those produced at 
RHIC. 

From a phenomenological perspective one of the most compelling evidences for the hydro- 
dynamic interpretation of RHIC flow results is the fact that the deviations from hydrody- 
namics are qualitatively reproduced by kinetic theory [HI, 13]. In particular, kinetic theory 
calculations generically reproduce the flattening of v^^pr) at higher p^, and the reduction of 
elliptic flow at large impact parameters. Some aspects of these kinetic theory results can be 
understood by considering the first viscous corrections to the thermal distribution function 



21[. These estimates motivated full viscous hydrodynamic simulations of the elliptic flow 
which will be performed in this work. Recently such viscous simulations have been performed 
by two other groups [22|, |23| and we will compare our results to these works in Section IVIBI 
A brief discussion of the history surrounding viscous relativistic hydrodynamics is given 
below. 
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B. Viscous Hydrodynamics 



The Navier-Stokes equations describe viscous corrections to ideal fluid flow by keeping 



terms up to first order in gradients of ideal quantities [24 
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The resulting equations are 
For instance, the stress tensor 
This is, of course, an unphysical 



parabolic which permit acausal signal propagation 
instantaneously adjusts to any thermodynamic force, diUj. 
picture since the stress tensor should relax to the thermodynamic forces over a typical 
collision timescale. 

One would therefore like a phenomenological theory that explains this relaxation correctly. 
Much work has been done in this direction but there is still no completely satisfactory theory. 
Probably the most used model is that of Israel and Stewart 26l . l27l| . but there are also others 



by Lindblom and Geroch [28|] , Pavon, Jou and Casas-Vasquez [291] and also by Ottinger and 
Grmela [H, [3o| which is used in this work. In fact a wide class of models was developed by 
Lindblom and Geroch in two separate papers 28j, |31 



All of the above theories have the same behavior: they relax on small time scales to 
the first-order relativistic Navier-Stokes equations and have some generalized entropy which 
increases as a function of time. It was shown by Lindblom 32] that for a large class class of 
these second order theories, the physical fields should be indistinguishable from the simple 
Navier-Stokes form. To paraphrase Lindblom; any measurement of the stress energy tensor 
or particle current on a time scale larger than the microscopic time scale will be indistin- 
guishable from the Navier-Stokes theory. The differences between the causal theories and 
the acausal Navier-Stokes equations are indicative of the corrections quantitatively captured 
by the full kinetic theory. Nevertheless, the causal theories provide a qualitative guide to 



the magnitude of these corrections [33|. However, the form of these corrections implicitly 



assumes a good quasi-particle description which may not exist in a strongly coupled plasma 
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There has been a large body of work in applying dissipative theories to central heavy-ion 
collisions [2l|, |35|, Perhaps a particle method will ultimately be the best way to include 
the effects of viscosity and the corresponding fluctuations in the stress tensor [3 71 . l38l . |39 



Even though the equations for non-central (2-|-l dimensions) dissipative hydrodynamics are 



known (e.g. 
collisions [il. 




only recently have results come out which simulate non-central heavy ion 
Further discussion of these results will be given in the discussion. 



II. THE HYDRODYNAMIC MODEL 



In the following section we outline the equations of motion for the hydrodynamical model 
used in the following simulations. We start by summarizing the well known first-order 
Navier-Stokes theory. Then we outline the equations required for a second-order causal 
description of dissipative fluid dynamics. This is done assuming a boost invariant expansion 
as first proposed by Bjorken [42|], where the equations of motion are expressed in terms of 
the proper time r = \/W—z^ and the spatial rapidity ?7s = |ln^. The cartesian coordinate 
z denotes the position along the beam axis while x,y label positions transverse to the beam 
axis. 
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A. 1^ Order Viscous Hydrodynamics - Navier Stokes 

Viscous hydrodynamics was originally outlined in the first-order Navier-Stokes approx- 
imation where the energy momentum tensor and baryon flux is a sum of their ideal and 
dissipative parts: 

T^"" = euV + {p + U)A'"' + n'"' , (2.1) 
n^' = nu^+j^, (2.2) 

where p, e, n and = (7, 7v) are the pressure, energy density, baryon density and four- 
velocity of the fluid. We use the convention that g'^'^ = diag(— 1, +1, +1, +1) and therefore 
u^u^ = — 1. The dissipative terms, vr and jd depend on the definition of the local rest frame 
(LRF) of the fluid. A specific form of tt^" and can be found using the Landau-Lifshitz 



definition [24] of the LRF (u^tt^'^ = 0), constraining the the entropy to increase with time 
and by working within the Navier-Stokes approximation (keeping terms to first order in 
gradients only) resulting in 

= _ry(V'^M'^ + V"m^ - ^A'^^VX) . (2-3) 
n = -CVpu^ , (2.4) 

fd = -(^)'V'^(f ) ' (2-5) 

where k, rj and ( are the heat conduction, shear and bulk viscosities of the fluid with tem- 
perature T and chemical potential /x. The viscous tensor is constructed with the differential 
operator = A^^^du where A^^" = g^^ + u^u^ is the local three-frame projector and 
d^u'' = d^u^ + r^ji^M'^ is the covariant derivative. 

The transport coefficients in a quark-gluon plasma and also in the hadronic gas were 
studied in Refs. It was found that the dominate dissipative mechanism was 



shear viscosity in both the QGP and hadronic gas. Bulk viscosity may however dominate 
in the transition region Heat transport can be ignored in the limit that ^ T which 
is the limit taken here. 

In the following work we will consider viscous effects in a quark-gluon plasma phase only. 
For this purpose we consider a constant shear to entropy ratio, rj/s = const and a massless 
gas p = l/3e. Future work will discuss viscosity in the mixed and hadronic phases. From 
this point on we will neglect the thermal conductivity. We keep the bulk viscosity in the 
equations for consistency, but always set C = in any calculations. 



B. 2 order Viscous Hydrodynamics 

In order to render a second order theory it is necessary to introduce additional variables. 
These variables will relax on very short time scales to the standard thermodynamic quantities 
in the first order theory, but an evolution equation for them is still required in order to avoid 
acausal signal propagation. One such theory that has been used in a number of works was 
introduced by Israel and Stewart j26j. Instead we use a theory developed by [H, Eoj due to 
its appealing structure when implemented numerically. However, as discussed above, all of 
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these theories should agree {i.e., they all relax on short time scales to the same the first-order 
equations) . 

We now summarize the evolution equations used in the current analysis following the 
mathematical structure outlined in Ref. 30(]. We use a simplified version of the model for 
deviations of the stress energy tensor close to equilibrium. The new dynamical variable that 
is introduced is the tensor variable c^u which will later be shown to be closely related to the 
velocity gradient tensor, 7?^,^. The tensor variable Cf^i, is conveniently defined to have the 
property 

c^^u" = , (2.6) 

and the energy momentum tensor is given by 

T'''^ = (e - u^f"'Pup)u^u'' + P^^ . (2.7) 

The explicit form of the stress tensor P^'^ is given in ^30( and has a fairly complicated 
form. The discussion in simplified by by considering small deviations from local thermal 
equilibrium and working in the local rest frame where the stress tensor is approximated as 

Ti^RF = PiS'' - , (2.8) 

where a is a small parameter related to the relaxation time (see appendix [C]) . The equations 
of motion are dictated by conservation of energy and momentum which is given by d^T^'^ = 0. 
In addition an evolution equation for the generalized mechanical force tensor is also needed 
and is given by ^3] 

u^{dxc^u - d^cxu - duC^x) = —Cf,u - —c^,u , (2.9) 

To T2 

where c and c are defined as the isotropic and traceless parts of the tensor variable c^t^ 
defined as 

Cf,u = ^{cl-l){v^,u + u^u^) , (2.10) 

C^U ~l~ U^Uy Cf^y -\- C^y . (2.1l) 

In the limit that the relaxation times (tq, T2) are very small the evolution equation yields 

d' = r2id,u' + dju' - '^5''dku'') + '^To5''dku' . (2.12) 

Substituting the above equation into Tf^jip and comparing the result to the Navier-Stokes 
equation (12. 5p the bulk and shear viscosities can be identified as 

rj = T2pa , 

C = lropa. (2.13) 



In the model proposed by Ottinger [30l| the quantity a is related to the equation of state, 
but in the linearized version it is simply treated as a constant parameter related to the 
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relaxation time. We fix a = 0.7 in all calculations, which then fixes the relaxation times 
{t2, tq) as a function of rj and (. The effects of varying a is shown in appendix I C 1[ 

It is natural to ask what is the effect of the relaxation time on the theory. In some sense 



this was already answered by Lindblom [32|. He showed that the physical fiuid must relax 
to a state that is indistinguishable from the Navier-Stokes form. Therefore we expect the 
physical velocity gradients to agree with those given by the auxiliary tensor variable c'^'^ as 
in Eq. ( ]2.12p . This is shown in Appendix ID 21 for various values of r]/s. We expect higher 
order gradient terms to be necessary when there are large deviations between any observable 
computed using the physical fields or the auxiliary field c^'^. This will be used as a gauge in 
order to find the limit of applicability of hydrodynamics 



1. 1+1 Dimensions 

We now outline the equations of motion for the stress-energy tensor and the generalized 
mechanical force tensor assuming a boost-invariant expansion as well as azimuthal symmetry 
with arbitrary transverse expansion. It is easiest to work in polar coordinates (r, r, 0, if) 
and since there is no dependence on or 77 the four-velocity can be expressed as = 
(7,7^^,0,0) where 7 = J- — In this coordinate system the metric tensor is given by, 

g^" = diag(-l, 1, l/r\ l/r^) 

The first two equations of motion are given by the conservation of energy and momentum, 
d^T^^ = for u = T and u = r. (Due to boost invariance and azimuthal symmetry the 
V = 7] and V = (j) equations are trivial.) 

9.T°° + 9,T°i = — (T°° + P33) _ lyoi ^2.14) 
r r 

a,r°^ -I- drT^^ = — - -iT^^ - p22) (2.15) 

r r 

where P^^ = r'^P'^^ and P^^ = t'^P^^. The evolution equations for the generic mechanical 
force tensor c'^'^ are: 

9,c" + vdrc^^ - -[(1 - c^^)dy + c^'^ay] = —c^^ - — (2.16) 

7 7x0 7x2 

o -22 I o -22 I -22 11 \ I 2 10 ""1-22 1 £22 /r, 
OrC + VOrC H (c — C ) + -C = C c (2-17) 

r r jTo 7x2 

drd'' + Vdrd'' + -(C33 + C°°) - -C^° = - —c'' (2.18) 

r r 7x0 7x2 

where c^^ = r^c^^ and c^^ = r^c^^. 



2. 1+2 Dimensions 

We now consider the 1-1-2 dimensional case without azimuthal symmetry but still having 
longitudinal boost invariance and use a coordinate system whereby the coordinates trans- 
verse to the beam axis are cartesian, {T,x,y,ri). Since there is no dependence on rj the 
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four-velocity can be expressed as = 7(1, Vx, Vy, 0) where 7 = , ^, In this coordinate 

system the metric tensor is given by, g'^'^ = diag(— 1, 1, 1, 

In this coordinate system the first three equations of motion are given by the u = t, x, 
and y components of the conservation law d^T^'^ = 0: 

9,T°° + 9,T°i + 9,T°2 = ^(r°o + r2p33) (2.19) 

drT^^ + 9^.r" + dyT^^ = ^T^° (2.20) 

S.T^o + dxT^' + dyT^' = ^T^o (2.21) 

The evolution equations for the generalized mechanical force tensor are: 



dyy' + 2[ic''-l)dxVx + c'H 

dyy+2[ic^^-l)dyVy + C^'c 

2 

12 I „12/a „, I Q „, \ I /„22 









7^2 


—22 


"22 
C^^ 




7^2 


— 1-33 


J_~33 




7T"2 




l)dyVx 






7^0 


7^2 



(2.22) 
(2.23) 
(2.24) 

(2.25) 



3. Initial Conditions 

The hydrodynamic simulation is a 2 + 1 boost invariant hydrodynamic model with an 
ideal gas equation of state p = | e. The temperature is related to the energy density with 
the Nf = 3 ideal QGP equation of state. We have chosen this extreme equation of state 
because the resulting radial and elliptic flow are too large relative to data on light hadron 
production. Thus, this equation of state will estimate the largest elliptic flow possible for a 
given shear viscosity. We note that for any non-central collision we have choosen a default 
impact parameter of b=6.5 fm. 

Aside from the equation of state, the hydrodynamic model is based upon reference [sf. 
At an initial time tq, the entropy is distributed in the transverse plane according to the 
distribution of participants for a Au-Au collision. Then one parameter, Cg, is adjusted to 
set the initial temperature and total particle yield. Specifically the initial entropy density 
in the transverse plane is 

s{x,y,To) = ^f^, (2.26) 
To ax ay 

where is the number of participants per unit area. The value Cs = 15 closely corresponds 
to the results of full hydrodynamic simulations [I, [g], 0] and corresponds to a maximum initial 
temperature of Tq = 420 MeV at impact parameter 6 = 0. With the entropy density specified 
the energy density can be determined. This requires inverting the equation of state. 
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In a viscous formulation we must also specify the viscous fields, i.e. the c'^'^ in the second 
order setup. Following the general philosophy outlined in Section [TO] we will choose the c^" 
such that the stress tensor deviations are 



fj,u 



-ri (V^u,) n = -CVX = (2.27) 



Since at time Tq the transverse flow velocity and the longitudinal flow velocity is Bjorken 
this means that at mid rapidity 

1 2 

VTxa; = TTyj; = - - 7T^^ = -r]dzU^ 11 = (2.28) 

To achieve this condition we first rewrite the flow equations for small c^i, and vanishing 
transverse flow. The Cij equations become 

drc'' = -_-_, (2.29) 

To T2 

a,c22 = -'—-'— ^ (2.30) 

To T2 
9 7^33 °33 

d^c^^ -t = -'—-'— . (2.31) 

In writing this we have used the fact that for small velocity ~ —u^vP . Then looking 
for the quasi stationary state we set the time derivatives to zero, and use the relations 
c*-' = \c\ 5^^ and = d^ + d^ to find that 

c" = ?^-|f (2.32) 



III. HYDRODYNAMIC RESULTS 



The equations outlined in the previous two sections were integrated numerically using the 
initial conditions described above. The algorithm and a discussion of the numerics can 
be found in appendix [Dl In this section we now show the results of the simulation. Before 
showing the results of the 2+1 dimensional simulation we outline some of the main physics 
points using results from the 1+1 dimensional case. 

Fig. [1] shows the energy density per unit rapidity (left) and the transverse velocity (right) 
at various times for both ideal hydrodynamics and for finite viscosity {rj/s = 0.2). The effect 
of viscosity is twofold. The longitudinal pressure is initially reduced and the viscous case 
does less longitudinal pdV work as in the simple Bjorken expansion [oj. This means that at 
early times the energy per rapidity decreases more slowly in the viscous case. The reduction 
of longitudinal pressure is accompanied by a larger transverse pressure. This causes the 
transverse velocity to grow more rapidly. The larger transverse velocity causes the energy 
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FIG. 1: (Color online) Plot of the energy density per unit rapidity (left) and of the transverse veloc- 
ity (right) at times of r = 1, 3, 6, 9 fm/c, for r]/ s = 0.2 (solid red line) and for ideal hydrodynamics 
(dotted blue line). 



density to deplete faster at later times in the viscous case. The net result is that a finite 
viscosity (even as large a.s rj/s = 0.2) does not integrate to give major deviations from the 
ideal equations of motion. A preliminary account of this effect was given long ago (46l |. 

We now present results of the 2+1 dimensional boost invariant hydrodynamic model. 
Fig. [2] shows contour plots of the energy density per unit rapidity in the transverse plane at 
proper times of r = 1, 3, 6, 9 fm/c. The initial conditions (r = 1) is taken from the Glauber 
model discussed before. 

Fig. [3] shows contour plots of the transverse velocity at the same times of r =1, 3, 6, 9 
fm/c. At r = 1 the figure is blank since the velocity in the transverse plane is zero as set by 
the initial conditions. By looking at the contours of constant v/c one can see that a finite 
viscosity increases the transverse velocity. 

Since we are interested in elliptic flow which originates from the initial spatial anisotropy 
of the collision region it is useful to see how the spatial and momentum anisotropy develop 



in time. We therefore look at the following three quantities |47 



(M) = ^;^y^ (3.1) 

where the double angular bracket ((■■•)) denote an energy density weighted average. The 
spatial ellipticity (e^.) is a measure of the spatial anisotropy as a function of time. The spatial 
anisotropy is what drives the momentum anisotropy (e^ ). T his quantity can be thought of as 



characterizing the p|, weighted integrated elliptic flow 48||. The final quantity {{vt)) is the 
average radial flow velocity. All three of these quantities are plotted in fig. Hlfor t]/s=0.2, 
0.05 and 10"^ 
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X (fm) X (fm) 

FIG. 2: (Color online) Contour plot of energy density per unit rapidity in the transverse plane. 
The contour values working outward are for r = 1 fm/c: 15, 10, 5, 1, 0.1, for r = 3 fm/c: 10, 5, 1, 
0.1, for T = 6 fm/c: 3, 2, 1, 0.5 and for r = 9 fm/c: 0.5, 0.375, 0.25, in units of GeV/fm^. 

As already shown in the 1+1 dimensional case the finite viscosity case does less longitu- 
dinal work. The longitudinal pressure is reduced while the transverse pressure is uniformly 
increased in the radial direction, i.e. gives no addition V2 component. This causes the 
transverse velocity (as seen in {{vt)) and fig. [3]) to grow more rapidly while lags behind 
the ideal case. Furthermore, the larger radial symmetric transverse velocity causes a faster 
decrease in the spatial anisotropy. This further frustrates the build-up of the momentum 
anisotropy e^. We therefore expect to see a decrease in the integrated V2 as the viscosity is 
increased. This is indeed the case as will be shown. However, this effect is small compared 
to the change in V2 from modifications of the off-equilibrium distribution function. 
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FIG. 3: (Color online) Contour plot of transverse velocity, v± = y + ^y- The inner most contour 
is for v± = 0.1 and increases in steps of Av± = 0.15. 



IV. FREEZEOUT 

As discussed in the introduction, ideal hydrodynamics is applicable when A^/p <^ L 
where L denotes the typical system size. When dissipative corrections are included, one 
must remember that the Navier Stokes equations are derived assuming that the relaxation 
time Tji is much smaller than the inverse expansion rate, r^d^u^ <^ 1. Therefore, in the 
simulations we determine the freezeout surface by monitoring the expansion rate relative to 
the relaxation time using a generalization of the freezeout criteria first proposed in [49|, 5^ 



and later in 51 
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FIG. 4: (Color online) Time evolution of the spatial ellipticity e^, the momentum anisotropy Cp, 
and the energy density weighted transverse flow {{v±)), see Eg. 13.11 
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Specifically, freezeout is signaled wlien^ 



- (4.1) 



'-d^u^' ~ TRdX ~ (4.2) 



This combination of parameters can be motivated from the kinetic theory estimates [52 . 
The pressure is p ~ e {v^-^) with {v^-^ the typical quasi-particle velocity and e the energy 
density. The viscosity is of order rj e (f^j^) with tr the relaxation time. Thus the 
freezeout condition is simply 

p / 

In the model we are considering rj/p = aT2 with a = 0.7 as described in Section fllBI 

The value of | can be considered as a parameter chosen to be smaller than one. The point 
is that as the above quantity becomes large the Navier Stokes approximation is no longer 
applicable and the simulation should freezeout. At this point one would need to include 
further higher order corrections in the gradients or switch to a kinetic approach. 

It is also convenient to have a definition for an analogous freezeout surface in the case of 
ideal hydrodynamics. One can think of keeping the freezeout surface fixed as rj/s is taken 
to zero. Dividing the freezeout criterion by rj/s and using s = {e + p)/T Ap/T we define: 

X = (4.3) 

which involves only quantities in the ideal simulation. This is a separate freezeout parameter 
independent of the viscosity. 

We show in fig. [5] contour plots of the freezeout surface for fixed x from both ideal (left 
plot) and viscous hydrodynamics (right plot). For fixed x the freezeout surfaces remain 
approximately the same in both cases. The freezeout surface from now on will be specified 
by X ill order to facilitate a comparison between the ideal and viscous cases when comparing 
spectra. 

We have typically chosen x ^^id rj/s m order that ^d^u^ = 0.6. Thus in Table [T] for 
rj/s = 0.2 we have x = 3.0 and ^d^u'^ = 0.6. However, for rj/s = 0.05 the freezeout 
parameter is x = 12 giving an unphysically large surface. This would normally not be the 
case in a more realistic model with a phase transition present, since in the hadronic phase 
the viscosity goes like rj ^ The change in scaling with temperature would cause the 
system to freezeout soon after hadronization. We plan on quantifying this statement in a 
future work. We therefore use {ri/p)dfj^u^ = 0.225 when rj/s = 0.05 giving x = 4.5. The thin 
solid curve in the right plot of fig. O shows this particular freezeout contour. In table [T] we 
summarize the freezeout parameters used throughout this work. For a given rj/s the most 
physical choice of freezeout parameter x is selected such that {ri/p)d^u^ ~ 0.6. However, 
if the viscosity becomes so small that the volume becomes unphysically large (such as for 
Tj/ s = 0.05) we set x = 4.5 as a maximum. These three physically motivated parameter sets 
are given in bold in the table. 

We should stress that the freezeout surface taken in this work is different from the typical 
constant temperature surface used in many hydrodynamic simulations. From fig. O one can 



In actual simulations we take {rj/p) dfj_u^ = 0.6 for most runs (see below). 



13 



see from the temperature map that the surface is not an isotherm and actually spans a 
very wide range of temperatures. The freezeout surface is understood by examining the 
expansion rate in Bjorken geometry 



u 

drU^ H h dxU^ + dyU^ . 

T 



The resulting surface is due to a competition between the first two terms in 
times and the last two terms at later times 



(4.4) 
at early 



77/s 




X 


0.05 


0.6 


12.0 


0.05 


0.225 


4.5 


0.05 


0.15 


3.0 


0.2 


0.9 


4.5 


0.2 


0.6 


3.0 


0.133 


0.6 


4.5 



TABLE I: Freezeout parameters used throughout this work. For a given rj/s the most physical 
choice of freezeout parameter x is selected such that {r}/p)d^u^ ~ 0.6. However, if the viscosity 
becomes so smah (such as for rj/s = 0.05) that the volume becomes unphysically large (see text for 
discussion) we set x = 4.5 as a maximum. These three physically motivated parameter sets are in 
bold. 
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FIG. 5: (Color online) Contour plot of various freezeout surfaces for central Au-Au collisions. Left: 
Surfaces from ideal hydrodynamics where the freezeout condition is set by the parameter x=l-5, 3 
and 4.5. Right: Corresponding viscous solution where r]/s was fixed by the condition ^d^u^ = 0.6. 
The thin solid black curve shows the contour set by ^d^u^ = 0.225 for comparison. 
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V. SPECTRA 



A. Anisotropy 

Before computing the differential spectrum we will compute the momentum anisotropy 
as a function of time. The momentum anisotropy A2 (which differs from f 2 by the placement 
of averages) is defined as 

^ {pD + {pI) Sn + S22' 

where S^^ is the sphericity tensor and can be related to the hydrodynamics fields (i.e. u'^ 
n) and moments of the ideal particle distribution function. The explicit form is given 
in appendix |X] and generalizes an appendix of Ollitrault to the viscous case. From a 
theoretical perspective, A2 is preferred because it is almost independent of the details of the 
particle content of the theory (ist . 

We plot A2 in the following manner. At a given proper time we integrate over the surface 
of constant which has developed by time r. The remaining part of the surface is fixed by 
integrating over the matter which has not frozen out (x < Xf.o.) at fixed proper time. This 
can be thought of as a freezeout surface with a fiat top at time r. As time moves forward 
eventually all of the matter is frozen out over a surface set by constant x yielding a constant 
A2. 

Figure [6] shows A2 for four different freezeout surfaces. The figure on the left shows 
the results using only the ideal contribution to the sphericity (regardless of if viscosity is 
present). This will be analogous to using only the ideal particle distribution function when 
generating the spectrum. First look at the solid black curves which are generated using 
ideal hydrodynamics and a specified x- Fo^' ^ larger value of x a larger space-time region is 
evolved by hydrodynamics producing a larger elliptic flow or A2. The true ideal case where 
hydrodynamics is universally applicable is given by x = 00. We see that for x = 4.5 most 
of the elliptic flow is reproduced. 

In order to assess the role of viscosity we first look at the figure on the left. The dashed 
curves show A2 for i]/s = 0.05 and Tj/s = 0.2 without including viscous corrections to the 
distribution function. (For clarity, these curves are shown only for x = 3.0 and x = 4.5.) 
Without the corrections to the distribution function the viscous corrections to A2 are modest. 
The right figure shows the analogous plot, this time including the viscous corrections to the 
distribution function. The corrections are much larger and we therefore expect the viscosity 
to decrease the integrated elhptic flow. 



B. Spectra 



The thermal jOr and differential V2 spectra of particles are generated using the Cooper- 
Frye formula 53|] given by 



E- 
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(5.2) 
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T (fm/c) T (fm/c) 

FIG. 6: (Color online) A2 (defined in Eq. I5.ip as a function of r. The solid black lines show the 
ideal result for x =1-5, 3.0, 4.5 and 00. Also shown in the right and left figures respectively are 
the viscous results with and without including the viscous correction to the distribution function, 
for X = 3.0 and 4.5 and ij/s = 0.2 (dashed green curve) and for t//s=0.05 (dotted blue curve). 



The thermal distribution function used in the Cooper-Frye formula above also needs to 
include corrections due to finite viscosity. We therefore write / = /o + 5/ where fo is 
the ideal particle distribution and 6f is the viscous correction which has been derived in 
appendix [B] and is given by 



Sf 



1 



2(e + p)T2 



/o(l + /o)p>^ 



TT 



(5.3) 



For boltzmann statistics /o(l + fo) is replaced by fo- The elliptic flow is defined as the 
weighted average of the yields with cos(20): 



v2{pt) = (cos(20)) 



J^^dcf) cos(2(/)) 



dN 



dypTdpTO 



Pt 



dN 



(5.4) 



dypTdpTd(f> 



where is the angle between the decaying particle's momentum (p^) and the azimuthal 
angle of the collision region. 

A typical freezeout surface for x = 3 at an impact parameter b=6.5 is shown in fig. [71 
Color gradients show the temperature profile on the freezeout surface and as noted before 
the surface is not necessarily an isotherm. 

Differential pt spectra for massless particles are shown in fig. [8] for two different freezeout 
surfaces: x = 3-0 (left) and x = 4.5 (right). In both plots the ideal case is shown by the 
solid red line. First we discuss changes to the spectra brought about by modifications to the 
equations of motion by looking at the spectra generated with the ideal particle distribution 
{fo only). For both values of viscosity and both freezeout choices a hardening of the spectra 
is observed. This is expected since viscosity tends to increase the transverse velocity. 

The effect from the viscous corrections to the distribution function are more subtle. At 
earlier times the transverse flow has not fully developed and the longitudinal pressure is 
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FIG. 7: (Color online) Freezeout surface for semi-central (b=6.5) Au-Au collisions for ry/s = 0.2 
and X = 3.0. 



reduced while the transverse pressure is increased |21|. This is a consequence of the fact 
that the shear tensor is traceless. The increase in transverse pressure leads to a hardening 
of the spectrum after integration over the space-time freezeout surface. This is the case 
for X = 3 even though the corrections are small. At later times the larger transverse flow 
alleviates some of the longitudinal shear. When the hydro is finally in a full 3D expansion, 
the viscous correction tends to reduce the transverse pressure. This changes the sign of the 
viscous correction term. This is seen for x = 4.5 where the viscous corrections soften the 
spectrum slightly. 

As discussed above, any observable created by using the auxiliary variable c^'^ should 
agree with the results using the physical velocity fields. Therefore we also show the viscous 
corrections calculated using the physical gradients (denoted by 5/g), i-e., in the local rest 
frame the ir^^ is approximated by 

7r*i = -r]{d'u^ + d^u' - , (5.5) 

3 

when computing 5f. 

Overall, the corrections to the spectra are small so it is hard to see any differences between 
the two calculations. This will not be the case for the differential elliptic flow where this 
comparison will be more important. 

Figure [9] shows the differential elliptic flow using the same parameter set from the pt 
spectrum. The solid red curves shows the ideal spectrum and, as expected, a larger elliptic 
flow is generated for x = 4.5 compared to x = 3 since a larger fraction of the space-time 
volume is described by hydrodynamics. 

The viscous correction to the equations of motion causes only a small change in the 
elliptic flow as seen by comparing the results at finite viscosity using fo only with the ideal 
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FIG. 8: (Color online) Differential transverse momentum spectra for Au-Au collisions at b=6.5 
fm. The left plot is for freeze-out parameter x = 3 and the right for x = 4.5. In both plots the 
ideal case is shown by the solid red curve. Then the viscous case is shown without including the 
viscous corrections to the distribution function and is denoted by fo- The addition of the viscous 
correction to the distribution function is generated in two different ways. Sf-,^ is calculated using 
the auxiliary tensor c^^ while Sfc is calculated using the physical gradients i.e., tt^'^ = —ri{d^d'^) . 



case. For x = 3 the change is almost negligible. For x 
2% at Pt = 2 GeV. 



4.5 deviations are on the order of 




FIG. 9: (Color online) Differential V2 spectra for Au-Au collisions at b=6.5 fm. The resulting 
curves are generated in the same way as described for the pT spectra in fig. [8] 

Including the viscous corrections to the distribution function can bring about large 
changes in the elliptic flow. We show the corrections due to the auxiliary variable by (J/tt 
and those from the gradients by Sfc and we expect the two results to agree. When the two 
results start to diverge the gradient expansion is no longer valid and a kinetic description is 
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really required. 

Based on our discussion in section [IV] the viscosity is what sets the freezeout surface. For 
ri/s = 0.2 we find that x = 3 (figure on left). In this case the viscous corrections are large 
but can only be trusted up to ~ 1 GeV. We also show for comparison the spectra for 
Tj/s = 0.05 which can be trusted past 2 GeV. For rj/s = 0.05 we take x = 4.5 for reasons 
discussed in section [IVl Again, the viscous correction decreases the elliptic fiow as a function 
of pt- Also shown are the spectra for rj/s = 0.2 and the corrections are larger. In both cases 
the spectra can be trusted past px = 2 GeV. 



VI. DISCUSSION AND COMPARISON WITH OTHER WORKS 



A. Discussion 



In summary we now make several conclusions regarding the effects of shear viscosity on 
heavy ion collisions. 

We first recall the setup. The paper is restricted to an ideal gas equation of state p = |e 
and sets the initial non-equilibrium fields to the value expected from the navier stokes 
equations vr*-' = —rj (d'^u^). The initial distribution of entropy density follows the distribution 
of p articipants. (This could be changed to a Color Glass Condensate model initial conditions 
.) The paper simulates a fluid model based on [30]] which is similar but differs from that 



of Israel and Stewart. However all models should ultimately agree on the magnitude of 
viscous corrections provided the viscosity is sufficiently small. 

Several technical notes warrant discussion here. An algorithm for a reliable solution of 
the viscous model was developed by Pareschi 45| and is presented in appendix [D] which 
achieves uniform numerical accuracy across a wide range of relaxation times. For small 
enough relaxation times the auxiliary fields vr*-' should relax to the form expected from 
the Navier-Stokes equation vr*-' ~ —i]{d'^u^). To see this good/reasonable convergence for 
small/modest viscosities see appendix ID 2[ Generically, relaxation models for viscosity have 
long time parameters (the shear viscosity rj in this case) and short time parameters. In the 
model considered here, a (see appendix [C]) is the short time parameter while in the Israel- 
Stewart theory this short time parameter is ti/[(€ + plr^r]. These short time parameters can 
be constrained by the /-sum rule |55|, l56|, l57| and is discussed further in appendix O In 
general the results should not depend on these short time parameters. 

We now summarize our physical results. The integrated viscous corrections to the fiow 
are small. This was seen in both the hydrodynamic fields and also in the differential and 
integrated elliptic fiow when the thermal distribution function was restricted to the ideal 
form. (The remainder of this paragraph discusses only results with this restriction.) For the 
integrated V2 this is seen in the left plot of fig. [6] where A2 is shown for ideal runs and viscous 
runs at 77/5 = 0.05 and 0.2. Corrections due to the modified fiow pattern are also small in 
the differential V2 spectrum as seen in fig. [9] by comparing the ideal and viscous runs (again 
with fo only.) Although there is the possibility for the elliptic fiow to be modified from 
variations in the freezeout surface across different runs this was minimized by freezing out 
on contours of constant x- O^is can see from fig. [5] that the space-time freezeout contours 
are about the same at zero and finite viscosity. The fact that only small changes in the fields 
are seen when including viscosity is not surprising. The time scale of any heavy ion collision 
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is much shorter then the time needed for dissipative effects to integrate and become large. 

Even though viscosity does not modify the flow strongly we have shown that there are 
still large corrections to the particle spectra due to off-equilibrium corrections to the ideal 
particle distribution function. Any bounds for the viscosity (at least from this paper) would 



have to come from the V2 spectra. As Lindblom [32j and earlier work by others [551] has 
clarified, any observable computed from the auxiliary fields vr*-' must agree with the same 
observable generated by the physical gradients —r]{d'^u^). When deviations are seen the 
viscous corrections can no longer be trusted. For a freezeout surface set by x = 4.5 the 
viscous corrections agree with gradients up to 2 GeV for viscosities as large as rj/s = 0.2 as 
seen in figure [H It is therefore safe to use only the auxiliary variable when generating spectra 
for this particular parameter set. In figure fTOl we show a summary plot of the differential 
elliptic flow. We now show one additional curve for rj/s = 0.133 yielding {ri/p)d^u^ = 0.6 
for this particular choice of freezeout surface. We believe that this choice of parameters is 
the closest physical scenario. The right plo t of figure [10] shows the measured elliptic flow 
as measured by the STAR collaboration [58|. We do not intend to make a comparison, but 
simply would like to keep the data in mind. Nevertheless since this simulation was performed 
with a massless gas which has the largest elliptic flow, it seems difficult to imagine that the 
r)/s > 0.35 will ever fit the data even if the initial conditions are modified along the lines of 
Ref. Q. 

Before a realistic comparison with data can be made the QGP/hadronic phase transition 
must be taken into account. In the vicinity of the phase transition it is possible that the shear 
viscosity may become very large. Also, a more realistic model for the hadronic gas would 
be the hard sphere model where 77 ~ ^. This would adjust at what point the simulation 
freezes out and would therefore effect spectrum. There is most likely a finite bulk viscosity 
due to the fluctuations of the QGP and hadron concentrations in the mixed phase or from 
chemical off-equilibrium in the hadronic phase 4J]. A final issue that should be taken into 
consideration is that particles of different mass could possibly freezeout on different surfaces. 
These issues will be addressed in a future work. 



B. Comparison 

We now compare our results to some other groups, first with the recent results of Song 



and Heinz [22| where they computed differential V2 spectrum in Cu-Cu collisions. 

One conclusion they found is that varying the initial conditions do not change the end 
result, even in the extreme condition when the equilibrium stress tensor is set to zero 
'^'^'^{'^0) = 0. This insensitivity is similar to the insensitivity to the short time parame- 
ter a indicated in appendix IC 1[ 

Song and Heinz also found that the viscosity substantially changes the flow. Their dif- 
ferential V2 spectra changes dramatically when viscosity is included, even if the particles 
freezeout using an ideal distribution. In our case there is almost no change in V2 when 
viscosity is present when freezing out with fo only, while in their case they see f 2 decrease 
by a factor of two at 3 GeV due to changes in the flow alone (see their fig. 4). It is possible 
that this difference is due to their inclusion of a phase transition in their equation of state. 
Once the plasma phase reaches the phase transition the momentum anisotropy (e^) stalls in 
their model staying constant for the entire mixed phase (and only growing slightly in the 
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FIG. 10: Left: Summary plot showing V2 for massless particles for simulations using ideal hydro 
and r]/s = 0.05,0.2. Right: Charged hadron V2 data using the standard reaction plane method as 
measured in Au-Au collisions at ^/s = 200 GeV for a centrality selection of 10% to 20% jssl- 



hadronic phase). Therefore after freezing out the elhptic flow is in some sense probing earher 
times where viscous corrections to the flow may be larger. Actually, if one looks at our ep in 
figure m the largest differences between the ideal and viscous cases is for earlier times. This 
explanation should be verified. 

When the viscous corrections to the distribution function are added we both see quali- 
tatively the same behavior. The viscous correction gets larger with transverse momentum 
eventually driving the elliptic flow below zero. It is impossible to make quantitative com- 
parisons until we include a phase transition and run simulations with smaller systems sizes, 
which we plan on doing in a future work. 



We now compare our work with that of Baier and Romatschke [59|, |60[ for the flow and 
Pt spectrum in central collisions. They also find that the viscosity does not integrate to 
significantly modify the ideal flow. We both flnd qualitatively the same behavior conflrming 



even earlier works [35|, |46|. Finite viscosity causes the temperature to drop slowly at earlier 
times and more quickly at later times compared to the ideal case. This effect was already 
discussed in section UTT] 

Although the freeze-out conditions in this work differ from those of Baier and Romatschke, 
we flnd qualitatively the same behavior when comparing pj- spectra, i.e. a hardening of 
spectra at large pr- 

In comparison with the differential V2 results by Romatschke and Romatschke 23| we 
see qualitatively the same behavior. In this case the comparison is more direct since both 
simulations were performed using Au-Au collisions at the same ^/s. Comparing our results 
for X = 4.5, we see that our V2 drops by ~50% at pt = 2 GeV. This is on the same order 



as seen in flg. 3 of Ref. [23|]. However they do not show contributions from flow and the 



distribution function separately so it is hard to make any deflnitive comparisons. 
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VII. CONCLUSIONS 



In this work we have outhned the equations of motion necessary for a casual description 
of viscous relativistic hydrodynamics and have shown results using initial conditions tuned 
to Au-Au collisions at y/s ~ 200 GeV. The results indicate that the viscous correction to 
the ideal equations of motion are small. The goal of this work was to calculate the viscous 
correction to differential V2 spectra. Even though modifications to V2 from the flow are 
small the effect of the off-equilibrium distribution function can bring about large changes 
in V2{pt)- By requiring observables calculated with the auxiliary fields to agree with those 
calculated with the physical gradients one can identify where a hydrodynamic description is 
reliable. 



APPENDIX A: SPHERICITY 

In this appendix we show (following jisj) how the momentum anisotropy, expressed 
through A2, can be related to hydrodynamic quantities. As discussed in the text A2 is 
defined as 

' 511 + ^22 {Pi) + {Piy ^ ^ 

The sphericity tensor Sij is calculated from the third moment of the momentum distri- 
bution function 

Sij = I da^S^-"", (A2) 

where da^ is a differential element of the freezeout surface and the third-rank tensor, S^^^ 
is defined as 

S^^P = I py'^pPf^p)-0-. (A3) 

In order to relate the sphericity tensor to hydrodynamic quantities we follow the same 
steps as was done in {4^ but also include the additional terms coming from viscous correc- 
tions. First one substitutes the expression for the momentum distribution function with the 
appropriate viscous correction term into the above equation for the third-rank sphericity 
tensor 

S^'^'P = + Si^""". (A4) 

The subscripts / and V correspond to the ideal and viscous contributions respectively and 
are defined as 

sr = /pw/o(p)^^^, (A5) 



'V 



2sT3 



2 



Sr'"^ (A6) 
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where we have defined the fifth-rank tensor 

Sr"^ = J pWp''p^fo{p)j0^- (A7) 

Lorentz invariance sets the form of both Sj and 6*5 as follows 

5f ^ = Au^u^uP + Big^^'uP + permutations) (A8) 
, 5^""^ = CuVu^u^'u^ + D{uVuPg''^ + permutations) + 

E{u^'g''Pg''^ + permutations). (A9) 



The coefficients A..E can be found in the same manner as was done previously in |48|. 
Quoting the result: 

A = n{2{{E'')) - m^) (AlO) 

B = n^ — ^ (All) 

C = n{lQ{{E^)) - lQm^{{E^)) + 3m^)/3 (A12) 
D = -n(8((E^)) - llm2((E2)) + 3m^)/15 (A13) 
E = n{{{E^)) - 2m^{{E^)) + m^)/15 (A14) 

The results for S}^ and S]"^ can be expressed in terms of hydrodynamic quantities: 

Sj^ = I {Aul + B)u>'d(r^ + 2B f u'da^ 

J (T J (7 

Sf = I {Aul + B)u^'da^ + 2B I u^day (A15) 
The result for Sy and can be found making use of the following three identities: 

U"'naf3 = 

^"^vr,/3 = 

u"Aa(3 = (A16) 
Therefore, the viscous correction to the sphericity tensor is given as 

n r 

-— / [SDulu^dcTf, + E {hu>'da^ + IQw'da^ - Qulu^'da^)] , (A17) 



and 



n f 

— / [3Dulu^da^ + E {hu^'da^ + Wday - Gm^m^c/ct^)] . (A18) 
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APPENDIX B: VISCOUS CORRECTION TO THE DISTRIBUTION FUNC- 
TION 



The thermal px and differential V2 spectra of particles are generated using the Cooper 
Frye formula |53|] 

J^N g 



E- 



11} 



Bl) 



where da^ is the normal vector to the freezeout surface set by the condition of constant x- 
For the geometry we are considering here we have 



p^u'^ = rriTU^ cosh(r]s) — piu^ — P2U^-, 
p^da^ = TirriT cos\i{r]s)daQ + p^dai + p'^da2)- 



(B2) 
(B3) 



Following |2ll . |6l[ we will make a second moment ansatz for the thermal distribution 
function. We first write the stress tensor as 



T^'^ = eu^u" + {p + n)A'"' + vr 



(B4) 



where tt'^^ is symmetric traceless and satisfies ix^'^u^ = 0. Then we subsequently make an 
ansatz for the thermal distribution f fo + 



5f 



1 



Ci C2 . 

— TT^;/ + — ff A^j, 



(B5) 



where C is a constant and the factor 1/ [(e + p)T^] has been inserted for later convenience. 
To determine the constant C we demand that 



d^p p^p^ 



(27r)3 

Working in the local rest frame this becomes a condition that 



d^p p^pip^p^- 



foil + fo) 



^1 I ^'^nx 



{e + p)T^ J (27r)3 

The integral over the three momentum in curly braces can be expressed as 
where 



d'^P I P I 



15(e + p)r2 J (27r)3 Ep 
Inserting Eq. (IBSp into Eq. (lB7p we see that 

1 



/o(l + /o). 



Ci — C2 



(B6) 

(B7) 

(B8) 
(B9) 

(BIO) 



We record two limiting cases of this integral. In the massless limit with zero chemical 
potential the integral is easily performed and yields 



90C(5) 



0.958. 
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(Bll) 



In the classical limit the factor /o(l + /o) is replaced by 

foil + fo) ^fo = e-(^p-^)/^, 



(B12) 



and the integral is easily performed using the integral representation of the modified Bessel 
functions 



1 



(e + p)T2 



irp3 



27r2 



1. 



(B13) 



APPENDIX C: RELAXATION TIME 



It was shown in [56|, |57[ that for a weakly interacting theory the transport time scale 
tr is much longer than the inverse temperature. This separation of time scales is seen in 
the spectral density which will have a sharp peak at small frequencies u; ~ ^ T. In 
[sgI . 13] a sum rule for this peak was derived relating the small time {t <^ tr) behavior of 
hydrodynamic correlators to the microscopic time scale. The statement of this sum rule can 
be written as: 



T 



a,x!,(k,t)UvA^T(e + p)(^) 



(CI) 



where v-p is p/i? and A is a cut-off such that 1 / tr <^ K <^ T . Xgg is the retarded correlator 
of and can be found in the framework of linear response theory. A small velocity field is 
turned on with a perturbing Hamiltonian 



c/W(x,t)T°^(x,t) 



(C2) 



and suddenly switched off at t=0: f*(x, t) = e^*6'(— t)t'Q(x). In the framework of linear 
response this yields 



9,(T°Xk,t)) = -x^Jk,tK(k). 



(C3) 



The stress tensor can be expressed as the equilibrium stress tensor plus small corrections: 



(T°°) = e + e(x,t) 
(T°*) = + (7Xx,t), 

where g = v(e + p). The linearized hydrodynamic equations are: 

dte + dig' = 
dtg^ + dir'^ = 0, 



(C4) 
(C5) 
(C6) 



(C7) 
(C8) 
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where 

r'^ = S'^p - n{d'u^ + d^u' - \5''diu^) - 5'Xdiu\ (C9) 

in the Navier-Stokes hmit. However, since we are interested in relating the short time 
parameters of the theory used in this work to microscopic quantities we take t^^ from the 
second-order equations 

t'^ =p{6'^ -ad^). (CIO) 

We now have all the pieces needed in order to evaluate the left hand side of eq. ICli First, 
for small c the evolution equation simplifies to 

dtc'^ - {d'u' + d'u') = -c\S'' + -id' - lc\5''). (Cll) 

To T2 6 

We can now differentiate eqn. lCSI with respect to time, substitute in the evolution equation 
and immediately take the t — > limit in order to obtain 

d^^g^\t=o = -c%dte - apd^pu' + d^u% (C12) 

We now make use of the first linearized hydrodynamic equation, dte = —{e+p)diV^, which 
can be substituted into eq. IC12[ After taking a spatial Fourier transform we get: 

d^g^t=o = -clie + p)kjkiv' - pal^v^ + k'k'v% (C13) 

This equation can be decomposed into its transverse and longitudinal pieces by defining 
gi = g^'T -\- ^g^ where k^g^^ = 0. Since we are only interested in the shear viscosity we can 
simply look at the transverse component 

d^g^^ = -pak\^, (C14) 
and when substituted into the sum rule, eqn. ICll we obtain the result 

a = ^. (C15) 





A similar analysis can be done for the Israel- Stewart equations [26|] as well: 

r^^ = pS'^ + ir'^ (C16) 
dtTi'^ = —l-riid'v^) - n'^] (C17) 



with the result 



r, = |. (C18) 
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FIG. 11: Differential V2 spectrum for Au-Au collisions at b=6.5 fm. Left: ij/s = 0.05 Right: 
rj/s = 0.2 Each figure shows spectrum calculated with the default small time parameter a = 0.7 
and half of this value, a = 0.35. 



1. Dependence on small time parameter 

It was discussed throughout the paper that the results should not depend on the small 
time parameter. In order to test this we have generated f 2 spectrum with a value of a = 0.35 
compared to the default value of a = 0.7 used throughout this work. This is shown in Fig.fTT] 
for rj/s = 0.5 (left) and rj/s = 0.2 (right) for a fixed freeze-out surface parameter % = 3. 
The ideal curves are also shown for reference. 

The deviations between the results using two different values of a are small. The flow 
is hardly changed as seen by comparing the spectrum generated using only the ideal dis- 
tribution function, /q. In this case, for both values of rj/s, the results differ by less than 
1%. When including the viscous correction to the spectrum the results still agree reasonably 
well. 



APPENDIX D: ALGORITHM 

In this appendix the algorithm used in order to solve Eqn's 12. 19112. 251 is outlined. The 
numerical evaluation of the above system of hyperbolic equations is difficult because one 
would like would like to achieve uniform numerical accuracy across a range of relaxation 
times. In order to achieve this we use a discretization method first proposed by Pareschi 



45l | which can numerically solve the above equations in both the stiff and unstiff regions. 



We use notation such that the term xf^ refers to the value of x at discrete time and 
grid point (x, y) = (xj, yj ) with i and j always referring to the x and y grid coordinates. Any 
variable absent of an index represents a continuous variable not yet specified at a given point. 
At times we use a simplified notation for the discretized fields such that u{xi,yj,t"') = ufj- 

For completeness we outline the integration routine developed by [il]. Our goal is to 
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solve equations of the form 



ut{x, y, t) + f^{u) + hy{u) = g{u), (Dl) 

where u, f, h, and g are arbitrary functions. We use the standard notation that for a finite 
volume element the value of a field at the point {x^^i,yj^i) and at a time t = t" is given 
by: 



<^m = ^j^^ u{x,y,ndxdy (D2) 

We then integrate eg. ID II over the region x [yj,yj_^i] x yielding: 



^"+1 = I 

Ax Ay 



/ / u{x,y,t'')dxdy + / {f{xi+i,y,t)-f{xi,y,t))dydt + 

J Xi Jvi J Vi J t" 



Xi yj yj 

Xi+i rxi+1 ryj+i /•t"+^ 

/ (^h{x,yj+i,t) — h{x,yj,t))dxdt + / / / g{x,y,t)dxdydt 

Xi Jf^ Jxi Jyj Jt" 



(D3) 



The first integral over the field u in the above equation fIDSP can be discretized by 
constructing a piecewise linear approximation of u(x,y,t) over the integration region: 



2 



u{x, y, r)dy = / [u{x, yj, t") + {y - yj)uy{x, yj, t")] dy + 
Vj Jyj 
rvj+i 

+ / Hx, yj+u r) + iy- yj+i)uyix, yj+i, t")] dy (D4) 



Using the corresponding linear approximation for the integration over [x, Xj+i] and per- 
forming the elementary integration over x and y the following discretization is found: 

7^ J J «(^' ^' = 4 i<j + "IVij + Kj+i + <+i,,+i) + 



Ax Ay 

For the time integrals over the fluxes (second and third term in eq. ID3I1 a general trape- 
zoidal rule is used: 



1 r^" 

— J^^ fix, y, t)dt ^ fifix, y, + uf{x, y, T) (D6) 
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where /"+° will be given explicitly by a predictor step to be defined later. The time 
integrals over the source term (last term in eq. ID3I will also be given by a general trapezoidal 
rule which will result in an implicit equation between the sources and charges. 



1 

At 



g{x, y, t)dt ^ Cgix, y, + V9ix, y, t") (D7) 



In order to ensure second order accuracy in space the second and third integrals over the 
fluxes in eq. ID3I are evaluated after the time integrations using the standard midpoint rule 
and trapezoidal rule depending on the time: 



for t = t""*"^ : J y, ~ , y, t) 



(D8) 



for t e [r, ■ ^ t)dx ^ i [h{xi+i, y, t) + h{xi, y, t)] (D9) 

After evaluating all the integrals in eq. ID3I using the above rules for discretization the 
final result for u^'^^ . , i is: 



«+2 J 



2 'J ^2 



At 



[//I f" _L f" 1 -L 7/1 f"+Q: _|_ fn+a rn+a rn+a \ _i_ 

t^yi,j^Ji,j+l Ji+l,j Ji+l,j+l) ^ '^\Ji,j ^ Ji,j+l Ji+lJ Ji+l,j+lJ\^ 



2A2/ 



J+2 2- 

(DIO) 



The terms at time t""*"" are taken from the solution of the predictor step: 



n+a _ n 



Ata(djl^^^ + dyhl^ - g^;^'^) (Dll) 



As shown in |45l], second order accuracy conditions give the weights used in the discretiza- 
tion as a function of a. We choose a = 1/3 with weights given by /i = —1/2, z/ = 3/2,^ = 
3/4, and 77 = 1/4. We can therefore rewrite the solution ID 101 in operator splitting form as: 
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I A^ [" ^ rn j_ rn j-ra rn \ , ( cn+a , cn+a rn+a rn+a W , 

I \n{ /?" + 7)'^ — /)" — /)" ^ + 7/^/)"+" + "^1 4- 

~ + {9ij + ^i+ij + 9i,j+i + + j+ij 

(D12) 

The solution of fl2.19H2.25l) are a coupled set of seven equations of the form IDll At each 
timestep the following steps are performed. Update the charges according to the first line 
of ID12I for u^^^ at each point on the grid. Then solve implicitly for u"'~^°' where the source 
terms are possibly functions of the additional six field equations. Next update the charges 
according to m*-^-' at each grid point. Do a final implicit solve for according to the 

final equation in lD12[ 

1. ID versus 2D 

In order to demonstrate the robustness of the above algorithm in two dimensions we com- 
pare the results from the 2D numerical solution for central collisions with the corresponding 
ID result. In fig. [T2]the solid red line shows the result of the energy density per unit rapidity 
from the ID case using rj/s = 0.2. For reference the ideal result is shown by the dotted blue 
line. The black points plotted on top of the red curve shows the corresponding result for 
the 2D case. The difference between the two cases is small as expected. The scatter in the 
black points gives a qualitative idea of the error due to the use of a rectangular grid. Fig. [12] 
shows the analogous figure for the transverse velocity. 

2. Gradients 

It was discussed in the text that Lindblom [32] found an important result regarding the 
form of the auxiliary tensor. In a large class of causal dissipative theories the physical fluid 
states must relax to a state that is indistinguishable from the Navier-Stokes form. The time 
scale that which this occurs in on the order of the microscopic particle interaction time. 

We therefore should check that the viscous stress tensor n'^'^ as computed from the aux- 
iliary tensor c'^" agrees with the stress tensor computed from the gradients of the veloc- 
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FIG. 12: Scatterplot of the energy density per unit rapidity (left) and of the transverse velocity 
(right). 

ity field. This is shown for various components of tt^'^ from simulations with viscosity of 
rj/s = 10-^ 0.05, 0.2 in figs. MUB 
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